########################################################################
##### This file replicates the first study (magazine) in the paper #####
########################################################################
rm(list = ls())
library(dplyr)
library(fixest)
library(stargazer)
library(texreg)

#########################################################################
### This code merges the propaganda and military data, preparing them ###
### for analysis in a ready-to-use format. Note that this block of ######  
### code may take considerable time to execute. Researchers intending ### 
### to proceed directly with the analyses can skip to line 110 to load ## 
### the resulting dataset from this section. ############################
#########################################################################

### load propaganda (magazine) data collected for this project
magazine <- read.csv("../datasets/0_magazine.csv", stringsAsFactors = FALSE)

### load violence data; calculate violence data by group and time
violence <- read.csv("../datasets/0_merged_violence.csv", stringsAsFactors = FALSE)

reformed <- data.frame()
allgrouptime <- unique(magazine[,c("group", "pretime3mon", "pretime2mon", "time")])
for (i in 1:nrow(allgrouptime)){
  grouptime <- allgrouptime[i,]
  # time span: pretime3mon to time
  temp <- violence[violence$group == grouptime$group &
                     (violence$date >= grouptime$pretime3mon & violence$date < grouptime$time),]
  totalevent3mon <- nrow(temp)
  statebased3mon <- sum(temp$attack_type == 1, na.rm = TRUE)
  totaldeaths3mon <- sum(temp$deaths, na.rm = TRUE)
  # time span: pretime2mon to time
  temp <- violence[violence$group == grouptime$group &
                     (violence$date >= grouptime$pretime2mon & violence$date < grouptime$time),]
  totalevent2mon <- nrow(temp)
  statebased2mon <- sum(temp$attack_type == 1, na.rm = TRUE)
  totaldeaths2mon <- sum(temp$deaths, na.rm = TRUE)
  civildeaths2mon <- sum(temp$deaths_civilians, na.rm = TRUE)

  reformed <- rbind(reformed,
                cbind(grouptime, totalevent3mon, statebased3mon, totaldeaths3mon,
                      totalevent2mon, statebased2mon, totaldeaths2mon, civildeaths2mon))
}

reformed$pretime3mon <- as.Date(reformed$pretime3mon)
reformed$pretime2mon <- as.Date(reformed$pretime2mon)
reformed$time <- as.Date(reformed$time)


## load violence data; reform violence data for the placebo test
reformedpost <- data.frame()
allgrouptime <- unique(magazine[,c("group", "posttime2mon", "time")])
for (i in 1:nrow(allgrouptime)){
  grouptime <- allgrouptime[i,]
  # time span: time to posttime2mon
  temp <- violence[violence$group == grouptime$group &
                     (violence$date > grouptime$time & violence$date <= grouptime$posttime2mon),]
  totalevent2monpost <- nrow(temp)
  statebased2monpost <- sum(temp$attack_type == 1, na.rm = TRUE)
  totaldeaths2monpost <- sum(temp$deaths, na.rm = TRUE)

  reformedpost <- rbind(reformedpost,
                cbind(grouptime, totalevent2monpost, statebased2monpost, totaldeaths2monpost))
}

reformedpost$posttime2mon <- as.Date(reformedpost$posttime2mon)
reformedpost$time <- as.Date(reformedpost$time)

magazine$time <- as.Date(magazine$time)
magazine$pretime3mon <- as.Date(magazine$pretime3mon)
magazine$pretime2mon <- as.Date(magazine$pretime2mon)
magazine$posttime2mon <- as.Date(magazine$posttime2mon)

# merge propaganda and violence
magazine <- merge(magazine, reformed, 
                  by = c("group", "pretime3mon", "pretime2mon", "time"), all.x = TRUE)
magazine <- merge(magazine, reformedpost, 
                  by = c("group", "posttime2mon", "time"), all.x = TRUE)

# calculate aggregated variables
magazine_agg <- magazine %>% 
  group_by(group, pretime3mon, pretime2mon, time, type, language, magazine, issue, women, rawtime, total_issue, year, competitors) %>%
  summarize(religscore = mean(religscore, na.rm = TRUE),
            posscore = mean(posscore, na.rm = TRUE),
            pagenum = max(as.numeric(gsub(".txt", "", page))))
magazine <- merge(magazine, magazine_agg[,c("group", "magazine", "issue", "language", "pagenum")], by = c("group", "magazine", "issue", "language"), all.x = TRUE)

### transform violent variables
magazine$logtotaldeaths3mon <- log(magazine$totaldeaths3mon + 1)
magazine$logstatebased3mon <- log(magazine$statebased3mon + 1)
magazine$logtotalevent3mon <- log(magazine$totalevent3mon + 1)
magazine$logtotaldeaths2mon <- log(magazine$totaldeaths2mon + 1)
magazine$logstatebased2mon <- log(magazine$statebased2mon + 1)
magazine$logtotalevent2mon <- log(magazine$totalevent2mon + 1)
magazine$logcivildeaths2mon <- log(magazine$civildeaths2mon + 1)

magazine$logtotaldeaths2monpost <- log(magazine$totaldeaths2monpost + 1)
magazine$logstatebased2monpost <- log(magazine$statebased2monpost + 1)
magazine$logtotalevent2monpost <- log(magazine$totalevent2monpost + 1)

write.csv(magazine, "../datasets/1_magazine.csv", row.names=FALSE)


########################################################################
########################## Regression Analyses #########################
########################################################################

magazine <- read.csv("../datasets/1_magazine.csv", stringsAsFactors = FALSE)

### Main Text Table 1: Military Power and Religiosity Score in Magazines
mod1 <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
mod1cl <- summary(mod1, cluster = ~issue)
mod1fe <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod1fecl <- summary(mod1fe, cluster = ~issue)

mod2 <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
mod2cl <- summary(mod2, cluster = ~issue)
mod2fe <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod2fecl <- summary(mod2fe, cluster = ~issue)

mod3 <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
mod3cl <- summary(mod3, cluster = ~issue)
mod3fe <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod3fecl <- summary(mod3fe, cluster = ~issue)

mod1noisis <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue + competitors | 
                      language + year, 
                    data = magazine[magazine$group != "Islamic State" &
                                      magazine$type %in% c("Magazine"),]) 
mod1clnoisis <- summary(mod1noisis, cluster = ~issue)
mod1fenoisis <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$group != "Islamic State" &
                                  magazine$type %in% c("Magazine"),]) 
mod1feclnoisis <- summary(mod1fenoisis, cluster = ~issue)

texreg(list(mod1cl, mod1fecl, mod2cl, mod2fecl, mod3cl, mod3fecl, mod1clnoisis, mod1feclnoisis), stars = c(0.05))


### Main Text Table 2: Placebo Test: Future Military Power and Past Religiosity Score in Magazines
mod1fe <- feols(religscore ~ logtotalevent2monpost + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod1fecl <- summary(mod1fe, cluster = ~issue)

mod2fe <- feols(religscore ~ logtotaldeaths2monpost + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod2fecl <- summary(mod2fe, cluster = ~issue)

mod3fe <- feols(religscore ~ logstatebased2monpost + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod3fecl <- summary(mod3fe, cluster = ~issue)

texreg(list(mod1fecl, mod2fecl, mod3fecl), stars = c(0.05))


### Main Text Table 3: Military Power and Tolerance of Secular Values
# Rescale the sentiment (positivity) score to make the coefficients more readable, as they would otherwise be too close to zero.
magazine$posscore2 <- magazine$posscore*100
# Define the secularity score as the inverse of the religiosity score to simplify interpretation.
magazine$secularscore <- -magazine$religscore

pos1 <- feols(posscore2 ~ logtotalevent2mon + women + frontpage + total_issue | 
                group + language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
pos1cl <- summary(pos1, cluster = ~issue)

pos2 <- feols(posscore2 ~ logtotalevent2mon + secularscore + women + frontpage + total_issue | 
                group + language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
pos2cl <- summary(pos2, cluster = ~issue)

pos3 <- feols(posscore2 ~ logtotalevent2mon*secularscore + women + frontpage + total_issue | 
                group + language + year, data = magazine[magazine$type %in% c("Magazine"),]) 
pos3cl <- summary(pos3, cluster = ~issue)

texreg(list(pos1cl, pos2cl, pos3cl), stars = c(0.05))


### Appendix Table A1: Summary of Magazines by Groups
descriptmag <- unique(magazine[,c("group", "time", "type", "language", "magazine_name", "magazine", "issue", "women", "year")])
descriptmag <- descriptmag[!(descriptmag$type %in% c("Newsletter")),]
unique(descriptmag$language)
descriptmag$language <- gsub("Arabic", "ar", descriptmag$language)
descriptmag$language <- gsub("Bengali", "bn", descriptmag$language)
descriptmag$language <- gsub("English", "en", descriptmag$language)
descriptmag$language <- gsub("French", "fr", descriptmag$language)
descriptmag$language <- gsub("German", "de", descriptmag$language)
descriptmag$language <- gsub("Indonesian", "id", descriptmag$language)
descriptmag$language <- gsub("Malay", "ms", descriptmag$language)
descriptmag$language <- gsub("Russian", "ru", descriptmag$language)
descriptmag$language <- gsub("Swahili", "sw", descriptmag$language)
descriptmag$language <- gsub("Turkish", "tr", descriptmag$language)
descriptmag$language <- gsub("Uighur", "ug", descriptmag$language)
descriptmag$language <- gsub("Urdu", "ur", descriptmag$language)
descriptmag$women <- ifelse(descriptmag$women == 1, "w", "m")

descriptmag2 <- descriptmag %>% 
  group_by(group, type, magazine, magazine_name, language, women) %>%
  summarize(from = min(as.numeric(year), na.rm = TRUE),
            to = max(as.numeric(year), na.rm = TRUE))

descriptmag2$all <- paste0(descriptmag2$magazine_name, " (", descriptmag2$from, "-", descriptmag2$to, ", ", descriptmag2$language, ", ", descriptmag2$women, ")")

descriptmag3 <- descriptmag2 %>% 
  group_by(group) %>%
  summarize(allmags = paste(all, collapse = ", "))

descriptmag3$group[descriptmag3$group == "Islamic State"] <- "Islamic State of Iraq and Syria"
descriptmag3$allmags <- gsub("\x90", "'", descriptmag3$allmags)

stargazer::stargazer(descriptmag3, summary = FALSE, rownames = FALSE)


### Appendix Table A4: Religiosity Score in Western and Non-Western Languages
westlang <- cbind(mean(magazine$religscore[magazine$type %in% c("Magazine") & magazine$westlang ==0]),
                  mean(magazine$religscore[magazine$type %in% c("Magazine") & magazine$westlang ==1]))
colnames(westlang) <- c("Non-Western Languages", "Western Languages")
stargazer(westlang, summary = FALSE)


### Appendix Table A6: Alternative Measures for Military Power
# calculate variable: number of publications
groupmagazine <- as.data.frame(table(magazine$group, magazine$magazine))
groupmagazine <- groupmagazine[groupmagazine$Freq!=0,]
groupmagazine <- cbind(names(table(groupmagazine$Var1)), table(groupmagazine$Var1))
groupmagazine <- unname(groupmagazine)
colnames(groupmagazine) <- c("group", "total_pubs")
magazine <- merge(magazine, groupmagazine, by = c("group"), all.x = TRUE)
magazine$total_pubs <- as.numeric(magazine$total_pubs)
magazine$total_pubs <- magazine$total_pubs * magazine$total_issue

# calculate variable: publication frequency
windows <- NULL
for (name in unique(magazine$magazine)){
  temp <- magazine[magazine$magazine==name,]
  window <- as.numeric(stringr::str_extract(diff(range(as.Date(temp$time)))/mean(temp$total_issue), "\\d+"))
  windows <- rbind(windows, c(name, window))
}
windows <- as.data.frame(windows)
colnames(windows) <- c("magazine", "avgtime")
windows$frequency <- ifelse(windows$avgtime==0, 0, 1/as.numeric(windows$avgtime))
magazine <- merge(magazine, windows, by = "magazine", all.x = TRUE)

# calculate variable: professionalism
colorfulness <- read.csv("../datasets/0_magazine_colorful.csv", stringsAsFactors = FALSE)
colorfulness$is_colorful <- as.numeric(colorfulness$is_colorful)
colnames(colorfulness)[colnames(colorfulness)=="subfolder"] <- "magazine"
colorfulness$language <- gsub("Magazine_", "", colorfulness$main_folder)
colorfulness$language[colorfulness$language=="Farsi_Pashto"] <- "Farsi"
colorfulness$language[!(colorfulness$language %in% magazine$language)]
magazine <- merge(magazine, colorfulness[colorfulness$magazine %in% magazine$magazine,2:4], 
                  by = c("magazine", "language"), all.x = TRUE)

mod_pubs <- feols(religscore ~ total_pubs + women + frontpage + total_issue | 
                    group + language + year, 
                  data = magazine[magazine$type %in% c("Magazine"),]) 
mod_pubscl <- summary(mod_pubs, cluster = ~issue)

mod_freq <- feols(religscore ~ frequency + women + frontpage + total_issue | 
                    group + language + year, 
                  data = magazine[magazine$type %in% c("Magazine"),]) 
mod_freqcl <- summary(mod_freq, cluster = ~issue)

mod_color <- feols(religscore ~ is_colorful + women + frontpage + total_issue | 
                     group + language + year, 
                   data = magazine[magazine$type %in% c("Magazine"),]) 
mod_colorcl <- summary(mod_color, cluster = ~issue)

texreg(list(mod_pubscl, mod_freqcl, mod_colorcl), stars = c(0.05))


### Appendix Table A7: Summary Statistics for Magazine Analyses
colnames(magazine)
sumstats <- rbind(summary(magazine$religscore[magazine$type %in% c("Magazine")]),
                  summary(magazine$posscore[magazine$type %in% c("Magazine")]),
                  summary(magazine$logtotalevent2mon[magazine$type %in% c("Magazine")]),
                  summary(magazine$logtotaldeaths2mon[magazine$type %in% c("Magazine")]),
                  summary(magazine$logstatebased2mon[magazine$type %in% c("Magazine")]),
                  summary(magazine$frontpage[magazine$type %in% c("Magazine")]),
                  summary(magazine$total_issue[magazine$type %in% c("Magazine")]),
                  summary(magazine$women[magazine$type %in% c("Magazine")]))
row.names(sumstats) <- c("Religiosity Score", "Sentiment Score", "Total Violent Events (log)", "Total Deaths (log)",
                         "State-based Battles (log)", "Front Page", "Total Issue", "Target Women")
stargazer(sumstats, summary = FALSE)


### Appendix Table A10: Military Power and Religiosity Score in Magazines (Without ISIS)
mod2noisis <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue + competitors | 
                      language + year, data = magazine[magazine$group != "Islamic State" &
                                                         magazine$type %in% c("Magazine"),]) 
mod2clnoisis <- summary(mod2noisis, cluster = ~issue)
mod2fenoisis <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                        group + language + year, 
                      data = magazine[magazine$group != "Islamic State" &
                                        magazine$type %in% c("Magazine"),]) 
mod2feclnoisis <- summary(mod2fenoisis, cluster = ~issue)

mod3noisis <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue + competitors | 
                      language + year, data = magazine[magazine$group != "Islamic State" &
                                                         magazine$type %in% c("Magazine"),]) 
mod3clnoisis <- summary(mod3noisis, cluster = ~issue)
mod3fenoisis <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$group != "Islamic State" &
                                  magazine$type %in% c("Magazine"),]) 
mod3feclnoisis <- summary(mod3fenoisis, cluster = ~issue)

texreg(list(mod2clnoisis, mod2feclnoisis, mod3clnoisis, mod3feclnoisis), stars = c(0.05))


### Appendix Table A11: Military Power and Religiosity Score in Magazines and Newsletters
mod1 <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod1cl <- summary(mod1, cluster = ~issue)
mod1fe <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod1fecl <- summary(mod1fe, cluster = ~issue)

mod2 <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod2cl <- summary(mod2, cluster = ~issue)
mod2fe <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod2fecl <- summary(mod2fe, cluster = ~issue)

mod3 <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue + competitors | 
                language + year, data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod3cl <- summary(mod3, cluster = ~issue)
mod3fe <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine", "Newsletter"),]) 
mod3fecl <- summary(mod3fe, cluster = ~issue)

mod1noisis <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue + competitors | 
                      language + year, 
                    data = magazine[magazine$group != "Islamic State" &
                                      magazine$type %in% c("Magazine", "Newsletter"),]) 
mod1clnoisis <- summary(mod1noisis, cluster = ~issue)
mod1fenoisis <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                        group + language + year, 
                      data = magazine[magazine$group != "Islamic State" &
                                        magazine$type %in% c("Magazine", "Newsletter"),]) 
mod1feclnoisis <- summary(mod1fenoisis, cluster = ~issue)

texreg(list(mod1cl, mod1fecl, mod2cl, mod2fecl, mod3cl, mod3fecl, mod1clnoisis, mod1feclnoisis), stars = c(0.05))


### Appendix Table A12: Military Power and Religiosity Score in Magazines (Three-month Window)
mod1 <- feols(religscore ~ logtotalevent3mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine"),]) 
mod1cl <- summary(mod1, cluster = ~issue)
mod1fe <- feols(religscore ~ logtotalevent3mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod1fecl <- summary(mod1fe, cluster = ~issue)

mod2 <- feols(religscore ~ logtotaldeaths3mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine"),]) 
mod2cl <- summary(mod2, cluster = ~issue)
mod2fe <- feols(religscore ~ logtotaldeaths3mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod2fecl <- summary(mod2fe, cluster = ~issue)

mod3 <- feols(religscore ~ logstatebased3mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine"),]) 
mod3cl <- summary(mod3, cluster = ~issue)
mod3fe <- feols(religscore ~ logstatebased3mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod3fecl <- summary(mod3fe, cluster = ~issue)

texreg(list(mod1cl, mod1fecl, mod2cl, mod2fecl, mod3cl, mod3fecl), stars = c(0.05))


### Appendix Table A13: Military Power and Religiosity Score in Magazines (Only Languages from the Non-Western World)
mod1 <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod1cl <- summary(mod1, cluster = ~issue)
mod1fe <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod1fecl <- summary(mod1fe, cluster = ~issue)

mod2 <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod2cl <- summary(mod2, cluster = ~issue)
mod2fe <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod2fecl <- summary(mod2fe, cluster = ~issue)

mod3 <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue + competitors | 
                language + year, 
              data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod3cl <- summary(mod3, cluster = ~issue)
mod3fe <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine") & magazine$westlang==0,]) 
mod3fecl <- summary(mod3fe, cluster = ~issue)

texreg(list(mod1cl, mod1fecl, mod2cl, mod2fecl, mod3cl, mod3fecl), stars = c(0.05))


### Appendix Table A14: Robust Standard Errors Clustered at the Group and Language Level

mod1fe <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod1fecl <- summary(mod1fe, cluster = ~group+language)

mod2fe <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod2fecl <- summary(mod2fe, cluster = ~group+language)

mod3fe <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$type %in% c("Magazine"),]) 
mod3fecl <- summary(mod3fe, cluster = ~group+language)

mod1fenoisis <- feols(religscore ~ logtotalevent2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$group != "Islamic State" & magazine$type %in% c("Magazine"),]) 
mod1feclnoisis <- summary(mod1fenoisis, cluster = ~group+language)

mod2fenoisis <- feols(religscore ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$group != "Islamic State" & magazine$type %in% c("Magazine"),]) 
mod2feclnoisis <- summary(mod2fenoisis, cluster = ~group+language)

mod3fenoisis <- feols(religscore ~ logstatebased2mon + women + frontpage + total_issue | 
                  group + language + year, 
                data = magazine[magazine$group != "Islamic State" & magazine$type %in% c("Magazine"),]) 
mod3feclnoisis <- summary(mod3fenoisis, cluster = ~group+language)

texreg(list(mod1fecl, mod2fecl, mod3fecl, mod1feclnoisis, mod2feclnoisis, mod3feclnoisis), stars = c(0.05))


### Appendix Table A19: Groups Demonstrate Triumph in Violent Actions; Do Not Express Remorse over Civilian Harm
mod1_power <- feols(power ~ logtotalevent2mon + women + frontpage + total_issue | 
                     group + language + year, 
                   data = magazine[magazine$type %in% c("Magazine"),]) 
mod1_powercl <- summary(mod1_power, cluster = ~issue)
mod2_power <- feols(power ~ logstatebased2mon + women + frontpage + total_issue | 
                     group + language + year, 
                   data = magazine[magazine$type %in% c("Magazine"),]) 
mod2_powercl <- summary(mod2_power, cluster = ~issue)
mod5_power <- feols(power ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                   group + language + year, 
                 data = magazine[magazine$type %in% c("Magazine"),]) 
mod5_powercl <- summary(mod5_power, cluster = ~issue)
mod6_power <- feols(power ~ logcivildeaths2mon + women + frontpage + total_issue | 
                      group + language + year, 
                    data = magazine[magazine$type %in% c("Magazine"),]) 
mod6_powercl <- summary(mod6_power, cluster = ~issue)

mod1_sad <- feols(sad ~ logtotalevent2mon + women + frontpage + total_issue | 
                      group + language + year, 
                    data = magazine[magazine$type %in% c("Magazine"),]) 
mod1_sadcl <- summary(mod1_sad, cluster = ~issue)
mod2_sad <- feols(sad ~ logstatebased2mon + women + frontpage + total_issue | 
                      group + language + year, 
                    data = magazine[magazine$type %in% c("Magazine"),]) 
mod2_sadcl <- summary(mod2_sad, cluster = ~issue)
mod5_sad <- feols(sad ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                      group + language + year, 
                    data = magazine[magazine$type %in% c("Magazine"),]) 
mod5_sadcl <- summary(mod5_sad, cluster = ~issue)
mod6_sad <- feols(sad ~ logcivildeaths2mon + women + frontpage + total_issue | 
                      group + language + year, 
                    data = magazine[magazine$type %in% c("Magazine"),]) 
mod6_sadcl <- summary(mod6_sad, cluster = ~issue)

texreg(list(mod1_powercl, mod5_powercl, mod2_powercl, mod6_powercl,
            mod1_sadcl, mod5_sadcl, mod2_sadcl, mod6_sadcl), stars = c(0.05))


### Appendix Table A20: Groups’ Openness in Discussing Death
mod_death <- feols(death ~ logtotaldeaths2mon + women + frontpage + total_issue | 
                    group + language + year, 
                  data = magazine[magazine$type %in% c("Magazine"),]) 
mod_deathcl <- summary(mod_death, cluster = ~issue)

mod_powerdeath <- feols(power ~ death + women + frontpage + total_issue | 
                        group + language + year, 
                      data = magazine[magazine$type %in% c("Magazine"),]) 
mod_powerdeathcl <- summary(mod_powerdeath, cluster = ~issue)

texreg(list(mod_deathcl, mod_powerdeathcl), stars = c(0.05))


### Appendix Section A6 Reverse Causality

## The section below restructures the data for Appendix Section A6 using 
## previously utilized page-level magazine data. It takes a considerable 
## amount of time to execute, so the resulting dataset has been saved to 
## ../datasets/1_miliprop.csv. Researchers can skip to line 610 to load 
## this dataset directly for replication. Uncomment the block below to 
## rerun the data preparation process.

# ## aggregate propaganda data at the group time level
# group_agg <- magazine[magazine$type %in% c("Magazine"),] %>% 
#   group_by(group, time) %>%
#   summarize(religscore = mean(religscore, na.rm = TRUE))
# group_agg$time <- as.Date(group_agg$time)
# group_agg$time_p1 <- group_agg$time + 30
# group_agg$time_p2 <- group_agg$time + 60
# group_agg$time_p3 <- group_agg$time + 90
# group_agg$time_p4 <- group_agg$time + 120
# group_agg$time_p5 <- group_agg$time + 150
# group_agg$time_p6 <- group_agg$time + 180
# group_agg$time_p7 <- group_agg$time + 210
# group_agg$time_p8 <- group_agg$time + 240
# 
# ### load violence data; calculate future military power for group and time at t+1, t+2, t+3
# violence <- read.csv("../datasets/0_merged_violence.csv", stringsAsFactors = FALSE)
# 
# miliprop <- data.frame()
# for (i in 1:nrow(group_agg)){
#   grouptime <- group_agg[i,]
#   # time span: t - t2
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time & violence$date < grouptime$time_p2),]
#   totalevent_t_t2 <- nrow(temp)
#   statebased_t_t2 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t_t2 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t1 - t3
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p1 & violence$date < grouptime$time_p3),]
#   totalevent_t1_t3 <- nrow(temp)
#   statebased_t1_t3 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t1_t3 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t2 - t4
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p2 & violence$date < grouptime$time_p4),]
#   totalevent_t2_t4 <- nrow(temp)
#   statebased_t2_t4 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t2_t4 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t3 - t5
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p3 & violence$date < grouptime$time_p5),]
#   totalevent_t3_t5 <- nrow(temp)
#   statebased_t3_t5 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t3_t5 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t4 - t6
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p4 & violence$date < grouptime$time_p6),]
#   totalevent_t4_t6 <- nrow(temp)
#   statebased_t4_t6 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t4_t6 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t5 - t7
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p5 & violence$date < grouptime$time_p7),]
#   totalevent_t5_t7 <- nrow(temp)
#   statebased_t5_t7 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t5_t7 <- sum(temp$deaths, na.rm = TRUE)
#   # time span: t6 - t8
#   temp <- violence[violence$group == grouptime$group &
#                      (violence$date >= grouptime$time_p6 & violence$date < grouptime$time_p8),]
#   totalevent_t6_t8 <- nrow(temp)
#   statebased_t6_t8 <- sum(temp$attack_type == 1, na.rm = TRUE)
#   totaldeaths_t6_t8 <- sum(temp$deaths, na.rm = TRUE)
#   
#   miliprop <- rbind(miliprop,
#                     cbind(grouptime,  totalevent_t_t2 = totalevent_t_t2,
#                           statebased_t_t2 = statebased_t_t2,
#                           totaldeaths_t_t2 = totaldeaths_t_t2,
#                           totalevent_t1_t3 = totalevent_t1_t3,
#                           statebased_t1_t3 = statebased_t1_t3,
#                           totaldeaths_t1_t3 = totaldeaths_t1_t3,
#                           totalevent_t2_t4 = totalevent_t2_t4,
#                           statebased_t2_t4 = statebased_t2_t4,
#                           totaldeaths_t2_t4 = totaldeaths_t2_t4,
#                           totalevent_t3_t5 = totalevent_t3_t5,
#                           statebased_t3_t5 = statebased_t3_t5,
#                           totaldeaths_t3_t5 = totaldeaths_t3_t5,
#                           totalevent_t4_t6 = totalevent_t4_t6,
#                           statebased_t4_t6 = statebased_t4_t6,
#                           totaldeaths_t4_t6 = totaldeaths_t4_t6,
#                           totalevent_t5_t7 = totalevent_t5_t7,
#                           statebased_t5_t7 = statebased_t5_t7,
#                           totaldeaths_t5_t7 = totaldeaths_t5_t7,
#                           totalevent_t6_t8 = totalevent_t6_t8,
#                           statebased_t6_t8 = statebased_t6_t8,
#                           totaldeaths_t6_t8 = totaldeaths_t6_t8))
# }
# 
# write.csv(miliprop, "../datasets/1_miliprop.csv", row.names = FALSE)

miliprop <- read.csv("../datasets/1_miliprop.csv", stringsAsFactors = FALSE)

miliprop$logtotalevent_t_t2 <- log(miliprop$totalevent_t_t2+1)
miliprop$logstatebased_t_t2 <- log(miliprop$statebased_t_t2+1)
miliprop$logtotaldeaths_t_t2 <- log(miliprop$totaldeaths_t_t2+1)
miliprop$logtotalevent_t2_t4 <- log(miliprop$totalevent_t2_t4+1)
miliprop$logstatebased_t2_t4 <- log(miliprop$statebased_t2_t4+1)
miliprop$logtotaldeaths_t2_t4 <- log(miliprop$totaldeaths_t2_t4+1)
miliprop$logtotalevent_t4_t6 <- log(miliprop$totalevent_t4_t6+1)
miliprop$logstatebased_t4_t6 <- log(miliprop$statebased_t4_t6+1)
miliprop$logtotaldeaths_t4_t6 <- log(miliprop$totaldeaths_t4_t6+1)
miliprop$logtotalevent_t6_t8 <- log(miliprop$totalevent_t6_t8+1)
miliprop$logstatebased_t6_t8 <- log(miliprop$statebased_t6_t8+1)
miliprop$logtotaldeaths_t6_t8 <- log(miliprop$totaldeaths_t6_t8+1)

### Appendix Table A9: Past Propaganda Does Not Significantly Affect Military Power

event_mod1 <- feols(logtotalevent_t_t2 ~ religscore | group, data = miliprop) 
event_mod2 <- feols(logtotalevent_t2_t4 ~ religscore | group, data = miliprop) 
event_mod3 <- feols(logtotalevent_t4_t6 ~ religscore | group, data = miliprop) 
event_mod4 <- feols(logtotalevent_t6_t8 ~ religscore | group, data = miliprop) 

deaths_mod1 <- feols(logtotaldeaths_t_t2 ~ religscore | group, data = miliprop) 
deaths_mod2 <- feols(logtotaldeaths_t2_t4 ~ religscore | group, data = miliprop) 
deaths_mod3 <- feols(logtotaldeaths_t4_t6 ~ religscore | group, data = miliprop) 
deaths_mod4 <- feols(logtotaldeaths_t6_t8 ~ religscore | group, data = miliprop) 

state_mod1 <- feols(logstatebased_t_t2 ~ religscore | group, data = miliprop) 
state_mod2 <- feols(logstatebased_t2_t4 ~ religscore | group, data = miliprop) 
state_mod3 <- feols(logstatebased_t4_t6 ~ religscore | group, data = miliprop) 
state_mod4 <- feols(logstatebased_t6_t8 ~ religscore | group, data = miliprop) 

texreg(list(event_mod1, deaths_mod1, state_mod1), stars = c(0.05))

### Appendix Figure A14: Past Propaganda Does Not Significantly Affect Military Power Across Different Time Windows
# extract coefficient and confidence intervals
results <- as.data.frame(cbind(c("event_mod1", "deaths_mod1", "state_mod1",
                   "event_mod2", "deaths_mod2", "state_mod2",
                   "event_mod3", "deaths_mod3", "state_mod3",
                   "event_mod4", "deaths_mod4", "state_mod4"),
                 sapply(list(event_mod1, deaths_mod1, state_mod1,
                             event_mod2, deaths_mod2, state_mod2,
                             event_mod3, deaths_mod3, state_mod3,
                             event_mod4, deaths_mod4, state_mod4), coef),
                 sapply(list(event_mod1, deaths_mod1, state_mod1,
                             event_mod2, deaths_mod2, state_mod2,
                             event_mod3, deaths_mod3, state_mod3,
                             event_mod4, deaths_mod4, state_mod4), se)))
colnames(results) <- c("model", "estimate", "std_error")
results$estimate <- as.numeric(results$estimate)
results$std_error <- as.numeric(results$std_error)
results$lower_ci <- results$estimate-1.96*results$std_error
results$upper_ci <- results$estimate+1.96*results$std_error

pdf("../outputplots/figureA14.pdf")
plot(NULL,
     main = NA,
     ylim = c(-0.1, 0.1), 
     xlim = c(0, 13), 
     ylab = "Effect of Propaganda on Future Military Power",
     xlab = "Interval: Propaganda to Military",
     cex.lab = 1,
     yaxt='n', xaxt = "n")
# Add axis labels for clarity
axis(1, at = c(2, 5, 8, 11), cex = 0.5, col.ticks = NA, 
     labels = c("0-2 months\nin the future", 
                "2-4 months\nin the future", 
                "4-6 months\nin the future", 
                "6-8 months\nin the future"))
axis(2, at = seq(-0.1, 0.1, by = 0.05), 
     col.ticks = 1, cex.axis = 1, las = 1)

# Add the null effect line
abline(h=0, col = "darkred")

# Plot each coefficient bar
colors <- c("purple", "lightcoral", "orange")
types <- c("Total Violent Events (log)", "Total Deaths (log)", "State-based Battles (log)")

for (i in 1:4) {
  points(i*3-1.5, results[i*3-2, 2], pch = 20, col = colors[1])
  segments(x0 = i*3-1.5,
           y0 = results[i*3-2, 4],
           x1 = i*3-1.5,
           y1 = results[i*3-2, 5], col = colors[1], lwd = 3)
  
  points(i*3-1, results[i*3-1, 2], pch = 20, col = colors[1])
  segments(x0 = i*3-1,
           y0 = results[i*3-1, 4],
           x1 = i*3-1,
           y1 = results[i*3-1, 5], col = colors[2], lwd = 3)
  
  points(i*3-0.5, results[i*3, 2], pch = 20, col = colors[1])
  segments(x0 = i*3-0.5,
           y0 = results[i*3, 4],
           x1 = i*3-0.5,
           y1 = results[i*3, 5], col = colors[3], lwd = 3)
}
# Add a legend
legend("topright", legend = types, pch = 16, col = colors, cex = 1)
dev.off()
